namespace NBody.Analysis
{
    using System;
    using DataFile;
    
    public class EnergyBinDMDEDJ : EnergyBin
    {
        Function x_axis;
        Function y_axis;
        Function slice;
        double slice_min, slice_max;
                
        public EnergyBinDMDEDJ()
        {
            n = 100;
            data = null;
            x_axis = new Function(E);
            y_axis = new Function(J2);
            slice = null;
        }
        
        public EnergyBinDMDEDJ(Function x, Function y, Function s, double min, double max)
        {
            n = 100;
            data = null;
            x_axis = x;
            y_axis = y;
            slice = s;
            slice_min = min;
            slice_max = max;
        }
        
        public override void Analyse(NBodySystem s)
        {
            Console.Error.WriteLine("Calculating dM/dEdJ");

            double xmax, xmin, ymax, ymin;
            calc_max(s, x_axis, out xmin, out xmax);
            calc_max(s, y_axis, out ymin, out ymax);
            xmax += 1e-10;
            ymax += 1e-10;
            
            Console.Error.WriteLine("  Xmin = {0:f6}", xmin);
            Console.Error.WriteLine("  Xmax = {0:f6}", xmax);
            Console.Error.WriteLine("  Ymin = {0:f6}", ymin);
            Console.Error.WriteLine("  Ymax = {0:f6}", ymax);
            
            double delta_x = (xmax - xmin) / (double)n;
            double delta_y = (ymax - ymin) / (double)n;
            //Console.Error.WriteLine("  DeltaX = {0}", delta_x);
            //Console.Error.WriteLine("  DeltaY = {0}", delta_y);
            
            data = new Table(n * n, 3, "Number density");
            
            for (int i = 0; i < n; i++)
            {
                for (int k = 0; k < n; k++)
                {
                    data[k + n * i, 0] = xmin + 0.5 * delta_x + i * delta_x;
                    data[k + n * i, 1] = ymin + 0.5 * delta_y + k * delta_y;
                    data[k + n * i, 2] = 0.0;
                }
            }
            
            for (int i = 0; i < s.NumParts; i++)
            {
                int x_pos = (int)((x_axis(s[i]) - xmin) / delta_x);
                int y_pos = (int)((y_axis(s[i]) - ymin) / delta_y);
                if (particle_in_slice(s[i]))
                    data[y_pos + n * x_pos, 2] += s[i].Mass / delta_x / delta_y;
            }
        }
        
        public static double E(Particle p)
        {
            return ((EnergyParticle)p).Energy;
        }
        
        public static double J2(Particle p)
        {
            return p.J2;
        }
        
        private void calc_max(NBodySystem s, Function f, out double min, out double max)
        {
            min = f(s[0]);
            max = f(s[0]);
            foreach (Particle p in s)
            {
                double val = f(p);
                if (val > max) max = val;
                if (val < min) min = val;
            }
        }
        
        private bool particle_in_slice(Particle p)
        {
            if (slice == null)
                return true;
            
            double val = slice(p);
            if (val < slice_max && val > slice_min)
                return true;
            else
                return false;
        }
    }
}
